GEO UNIV’R Tunisie 2024

LUN3 - Acquisition de données géographiques et visualisations de base

Auteur·rice

Ronan Ysebaert, Nicolas Lambert, Elina Marveaux

Date de publication

2024-05-13

Note

Ce support est inspiré en de nombreux points du manuel Géomatique avec R (Giraud et Pecout 2024). Et on les en remercie chaleureusement 🙏.

Ce document montre les solutions techniques pour importer, explorer, manipuler, visualiser (simplement) et exporter des données géographiques avec R :

Les exemples proposés sont tantôt adaptés aux contexte tunisien ou au contexte africain.

Packages utilisés

  • readxl : importer des fichiers Excel.
  • plotly : créer des graphiques interactifs.
  • sf : importer, manipuler et exporter des données géographiques vectorielles.
  • terra : importer, manipuler et exporter des fichiers raster.
  • leaflet : réaliser une carte interactive.
  • geodata : accéder à des jeux de données géographiques de référence dans le monde.
  • rnaturalearth : accéder à des fonds de carte du Monde.
  • wbstats : utiliser les jeux de données de la Banque Mondiale via son API.
  • osmextract : télécharger des données OpenStreetMap.
  • osrm : calculer des temps de parcours routiers via l’engin de routage OSRM.

Si vous n’avez pas ces packages, vous pouvez les installer en exécutant la ligne suivante dans la console.

install.packages(c('readxl',  'sf', 'terra', 'leaflet', 'geodata',  'wbstats', 
                   'rnaturalearth', 'osmextract', 'osrm', 'plotly'))


From scratch

En partant de rien et avec quelques lignes de code on peut créer un jeu de données spatial pour en faire, par exemple, une carte interactive.

Trois vecteurs sont créés, puis transformés en data.frame et enfin convertis en objet spatial avec la fonction st_as_sf() du package sf. Cela permet de créer rapidement une carte interactive avec le package leaflet.

# 3 vecteurs avec nom, longitude et latitude
name <- c("Sentido Bellevue Park", "Université de Sfax")
long <- c(10.579, 10.742)
lat <- c(35.913, 34.736)

# Transformer en data.frame
loc <- data.frame(name, long, lat)

# Transformer en objet spatial
library(sf)
loc <- st_as_sf(loc, 
                coords = c("long", "lat"),
                crs = 4326)

# Visualiser simplement avec leaflet
library(leaflet)

leaflet(loc) |> 
  addTiles() |> 
  addCircleMarkers(popup = loc$name) # Hôpitaux


Tableaux de données

Fichiers .csv

Import des données des délégations au format csv. On lui donne le nom de del_df.

del_df <- read.csv("data/tun/don_del.csv", sep = ";", dec = ",")

Avec R, il est très important de maîtriser la nature des objets importés ou transformés, et les convertir le cas échéant. La fonction class() nous permet de constater qu’il s’agit d’un data.frame. La fonction str() permet de détailler son contenu.

class(del_df)
[1] "data.frame"
str(del_df)
'data.frame':   264 obs. of  30 variables:
 $ del_code    : chr  "TN.AN.AR" "TN.AN.ET" "TN.AN.KA" "TN.AN.LS" ...
 $ del_nom_fr  : chr  "Ariana Ville" "Cité Ettathamen" "kalaât El Andalous" "Soukra" ...
 $ del_nom_ar  : chr  "أريانة المدينة" "حي التضامن" "قلعة الاندلس" "سكرة" ...
 $ gou_code    : chr  "AN" "AN" "AN" "AN" ...
 $ gou_nom     : chr  "Ariana" "Ariana" "Ariana" "Ariana" ...
 $ gou_cap     : int  1 0 0 0 0 0 0 1 0 0 ...
 $ gou_cap_dist: num  1 8.8 19.8 6.7 8 8.8 15 1 9.5 7.4 ...
 $ reg_code    : chr  "NE" "NE" "NE" "NE" ...
 $ reg_nom     : chr  "Nord-est" "Nord-est" "Nord-est" "Nord-est" ...
 $ popto_2004  : int  97687 78311 23045 89151 53752 60896 19404 32329 27977 31792 ...
 $ popco_2004  : int  97687 78311 15313 89151 40176 53911 8909 32329 27977 31792 ...
 $ immig_2004  : int  16961 5651 728 19129 7258 14053 1298 3963 5831 5508 ...
 $ emigr_2004  : int  15426 5245 528 2832 985 1443 723 9381 1135 3991 ...
 $ mobil_2004  : int  32387 10896 1256 21961 8243 15496 2021 13344 6966 9499 ...
 $ menag_2004  : int  27468 11950 4709 21590 17119 14276 4215 7941 6498 8462 ...
 $ ordin_2004  : int  9751 430 188 2979 685 2327 181 834 1176 1963 ...
 $ porta_2004  : int  22524 4505 1865 13321 7087 9022 1775 4908 4178 6397 ...
 $ telef_2004  : int  18596 3824 829 9262 7019 6339 1176 4002 2755 4747 ...
 $ popto_2014  : int  114486 84312 26796 129693 89884 106414 24503 31128 40101 34962 ...
 $ popco_2014  : int  114486 84312 18211 129693 58641 94961 11351 31128 40101 34962 ...
 $ immig_2014  : int  15637 5028 1104 20786 14400 20128 1997 2392 5855 3957 ...
 $ emigr_2014  : int  20448 6752 752 5528 1828 2521 822 8495 1883 3803 ...
 $ mobil2014   : int  36085 11780 1856 26314 16228 22649 2819 10887 7738 7760 ...
 $ menag_2014  : int  32498 22087 6554 33981 22781 27574 5922 8506 10066 9996 ...
 $ ordin_2014  : int  25474 6836 1701 18191 8326 14284 1723 4339 5494 6471 ...
 $ porta_2014  : int  32308 21715 6305 33683 22531 27338 5811 8363 9970 9901 ...
 $ telef_2014  : int  19942 3496 428 9831 4249 7505 637 2918 3110 4901 ...
 $ popto_2010  : num  109500 82922 24367 107829 69247 ...
 $ surfa_2010  : num  18.61 3.38 188.21 27.89 22.91 ...
 $ idr_2011    : num  0.638 0.386 0.383 0.557 0.466 0.531 0.358 0.489 0.51 0.523 ...

On visualise les premières lignes avec la fonction head().

head(del_df, 5)
  del_code         del_nom_fr     del_nom_ar gou_code gou_nom gou_cap
1 TN.AN.AR       Ariana Ville أريانة المدينة       AN  Ariana       1
2 TN.AN.ET    Cité Ettathamen     حي التضامن       AN  Ariana       0
3 TN.AN.KA kalaât El Andalous   قلعة الاندلس       AN  Ariana       0
4 TN.AN.LS             Soukra           سكرة       AN  Ariana       0
5 TN.AN.MN             Mnihla       المنيهلة       AN  Ariana       0
  gou_cap_dist reg_code  reg_nom popto_2004 popco_2004 immig_2004 emigr_2004
1          1.0       NE Nord-est      97687      97687      16961      15426
2          8.8       NE Nord-est      78311      78311       5651       5245
3         19.8       NE Nord-est      23045      15313        728        528
4          6.7       NE Nord-est      89151      89151      19129       2832
5          8.0       NE Nord-est      53752      40176       7258        985
  mobil_2004 menag_2004 ordin_2004 porta_2004 telef_2004 popto_2014 popco_2014
1      32387      27468       9751      22524      18596     114486     114486
2      10896      11950        430       4505       3824      84312      84312
3       1256       4709        188       1865        829      26796      18211
4      21961      21590       2979      13321       9262     129693     129693
5       8243      17119        685       7087       7019      89884      58641
  immig_2014 emigr_2014 mobil2014 menag_2014 ordin_2014 porta_2014 telef_2014
1      15637      20448     36085      32498      25474      32308      19942
2       5028       6752     11780      22087       6836      21715       3496
3       1104        752      1856       6554       1701       6305        428
4      20786       5528     26314      33981      18191      33683       9831
5      14400       1828     16228      22781       8326      22531       4249
  popto_2010 surfa_2010 idr_2011
1     109500     18.612    0.638
2      82922      3.376    0.386
3      24367    188.206    0.383
4     107829     27.895    0.557
5      69247     22.907    0.466

Fichiers Excel

La package readxl permet l’import de fichiers Excel. La fonction read_excel a plusieurs arguments utiles pour spécifier en entrée le format des colonnes (texte, numérique), ne pas considérer les n premières lignes du fichier, etc.

Nous importons ici un fichier Excel dérivé du World Population Prospects des Nations Unies, qui met à disposition toute sorte d’indicateurs démographiques (population, natalité, mortalité). On dispose ici des dénombrements de population par tranche d’âge quinquennale de 1950 à 2021 (feuille estimates).

library(readxl)

world <- read_excel("data/world/WPP2022_POP_F02_1_POPULATION_5-YEAR_AGE_GROUPS_BOTH_SEXES.xlsx", 
                 sheet = "Estimates",
                 skip = 16,
                 col_types = c(rep("text", 10), rep("numeric", 22))) 

Manipulation de data.frame

Dans les faits ces fichiers nécessitent généralement d’être reformatés et réorganisés pour pouvoir être interprétés dans un logiciel quel qu’il soit… Commençons par comprendre la structure des données : noms de colonnes, années disponibles etc.

# Conversion en data.frame
world <- data.frame(world)

# Nom des colonnes
colnames(world)
 [1] "Index"                               
 [2] "Variant"                             
 [3] "Region..subregion..country.or.area.."
 [4] "Notes"                               
 [5] "Location.code"                       
 [6] "ISO3.Alpha.code"                     
 [7] "ISO2.Alpha.code"                     
 [8] "SDMX.code.."                         
 [9] "Type"                                
[10] "Parent.code"                         
[11] "Year"                                
[12] "X0.4"                                
[13] "X5.9"                                
[14] "X10.14"                              
[15] "X15.19"                              
[16] "X20.24"                              
[17] "X25.29"                              
[18] "X30.34"                              
[19] "X35.39"                              
[20] "X40.44"                              
[21] "X45.49"                              
[22] "X50.54"                              
[23] "X55.59"                              
[24] "X60.64"                              
[25] "X65.69"                              
[26] "X70.74"                              
[27] "X75.79"                              
[28] "X80.84"                              
[29] "X85.89"                              
[30] "X90.94"                              
[31] "X95.99"                              
[32] "X100."                               
# Années disponibles
unique(world$Year)
 [1] 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964
[16] 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979
[31] 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994
[46] 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009
[61] 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021   NA

On peut alors filtrer les dimensions du tableau, calculer des indicateurs, etc. qui pourront être utilisés pour des statistiques de base, des graphiques voire de futures représentations cartographiques.

# Somme de plusieurs colonnes
world$POP_TOT <- rowSums(world[,c(12:32)]) # Population totale

# Renommer une colonne
names(world)[3] <- "Name"

# Créer des indicateurs
# Part des jeunes
world$YOUNG_RT <- (world$X0.4 + world$X5.9 + world$X10.14) / world$POP_TOT * 100 
# Part des personnes âgées
world$OLD_RT <- (world$X65.69 + world$X70.74 + world$X75.79 + world$X80.84 +
                   world$X85.89 + world$X90.94 + world$X95.99) / world$POP_TOT * 100 

# Sélection d'une ligne
df <- world[world$Year == 2021,] # Année de référence

# Sélectionner plusieurs lignes
afr <- c("910", "911", "912", "913", "914")
df <- df[df$Parent.code %in% afr,] # Pays africains

# Ordonner selon les valeurs d'une colonne
df <- df[order(df$YOUNG_RT, decreasing = TRUE),]

Il existe de nombreuses fonctions pour l’analyse statistique avec R. La plus basique étant probablement summary().

# Résumé stat
summary(df)
    Index             Variant              Name              Notes          
 Length:58          Length:58          Length:58          Length:58         
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
 Location.code      ISO3.Alpha.code    ISO2.Alpha.code    SDMX.code..       
 Length:58          Length:58          Length:58          Length:58         
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
     Type           Parent.code             Year           X0.4         
 Length:58          Length:58          Min.   :2021   Min.   :    0.22  
 Class :character   Class :character   1st Qu.:2021   1st Qu.:  299.04  
 Mode  :character   Mode  :character   Median :2021   Median : 1971.23  
                                       Mean   :2021   Mean   : 3576.71  
                                       3rd Qu.:2021   3rd Qu.: 4280.10  
                                       Max.   :2021   Max.   :34830.80  
      X5.9               X10.14             X15.19              X20.24         
 Min.   :    0.275   Min.   :    0.26   Min.   :    0.197   Min.   :    0.166  
 1st Qu.:  287.524   1st Qu.:  260.32   1st Qu.:  231.017   1st Qu.:  218.739  
 Median : 1799.815   Median : 1596.25   Median : 1340.744   Median : 1073.518  
 Mean   : 3237.939   Mean   : 2880.72   Mean   : 2481.160   Mean   : 2123.874  
 3rd Qu.: 3812.367   3rd Qu.: 3456.02   3rd Qu.: 3054.270   3rd Qu.: 2727.918  
 Max.   :30567.689   Max.   :26974.49   Max.   :22929.125   Max.   :18806.200  
     X25.29              X30.34              X35.39              X40.44        
 Min.   :    0.252   Min.   :    0.244   Min.   :    0.287   Min.   :   0.267  
 1st Qu.:  201.104   1st Qu.:  188.732   1st Qu.:  152.122   1st Qu.: 122.764  
 Median :  960.347   Median :  880.112   Median :  738.572   Median : 573.226  
 Mean   : 1858.879   Mean   : 1643.911   Mean   : 1421.126   Mean   :1171.518  
 3rd Qu.: 2325.300   3rd Qu.: 1923.398   3rd Qu.: 1581.452   3rd Qu.:1374.119  
 Max.   :15675.790   Max.   :13171.290   Max.   :11591.915   Max.   :9893.223  
     X45.49             X50.54             X55.59             X60.64        
 Min.   :   0.358   Min.   :   0.486   Min.   :   0.452   Min.   :   0.496  
 1st Qu.:  99.267   1st Qu.:  84.229   1st Qu.:  67.378   1st Qu.:  54.439  
 Median : 469.567   Median : 374.635   Median : 291.823   Median : 233.514  
 Mean   : 938.585   Mean   : 770.335   Mean   : 616.303   Mean   : 473.997  
 3rd Qu.:1124.716   3rd Qu.: 851.100   3rd Qu.: 667.341   3rd Qu.: 534.762  
 Max.   :7846.818   Max.   :6096.186   Max.   :4876.240   Max.   :3778.750  
     X65.69             X70.74              X75.79              X80.84       
 Min.   :   0.399   Min.   :   0.4345   Min.   :   0.3185   Min.   :  0.167  
 1st Qu.:  40.704   1st Qu.:  26.6259   1st Qu.:  16.2689   1st Qu.:  9.041  
 Median : 156.874   Median :  94.0440   Median :  58.9380   Median : 29.567  
 Mean   : 349.206   Mean   : 235.3130   Mean   : 140.6908   Mean   : 70.866  
 3rd Qu.: 377.993   3rd Qu.: 237.7533   3rd Qu.: 132.5434   3rd Qu.: 68.634  
 Max.   :2757.823   Max.   :1819.8435   Max.   :1102.5085   Max.   :494.341  
     X85.89            X90.94           X95.99            X100.         
 Min.   :  0.079   Min.   : 0.038   Min.   : 0.0075   Min.   :0.000000  
 1st Qu.:  4.053   1st Qu.: 0.968   1st Qu.: 0.0685   1st Qu.:0.004375  
 Median :  9.991   Median : 2.052   Median : 0.3372   Median :0.031750  
 Mean   : 28.099   Mean   : 7.973   Mean   : 1.4798   Mean   :0.219655  
 3rd Qu.: 27.383   3rd Qu.: 7.602   3rd Qu.: 1.2139   3rd Qu.:0.152125  
 Max.   :178.652   Max.   :81.603   Max.   :22.9160   Max.   :3.797500  
    POP_TOT            YOUNG_RT         OLD_RT      
 Min.   :     5.4   Min.   :13.97   Min.   : 1.677  
 1st Qu.:  2388.4   1st Qu.:34.33   1st Qu.: 2.696  
 Median : 12774.0   Median :40.41   Median : 3.152  
 Mean   : 24028.9   Mean   :38.12   Mean   : 4.298  
 3rd Qu.: 28556.3   3rd Qu.:43.45   3rd Qu.: 4.161  
 Max.   :213401.3   Max.   :48.90   Max.   :26.714  

Représentations graphiques

Les possibilités offertes en terme de représentations graphiques sont nombreuses avec R ! Car c’est aussi un langage qui a été créé pour cela.

De base

En une ligne de code on peut créer des représentations graphiques variées, adaptées à la nature statistique des variables à explorer.

par(mar = c(2,2,2,2), mfrow = c(2, 3))

barplot(df$YOUNG_RT, main = "Diagramme en barre")
boxplot(df$YOUNG_RT, main = "Boîtes à moustache")
hist(df$YOUNG_RT, main = "Histogrammes")
hist(df$YOUNG_RT, freq = FALSE, main = "Histogrammes et densité")
lines(density(df$YOUNG_RT), col = "blue")
stripchart(df$YOUNG_RT, method = "jitter", pch = 16,
           main = "Diagrammes de dispersion")
plot(data = df, YOUNG_RT ~ OLD_RT, main = "Nuages de points")

Ces graphiques sont paramétrables avec une série d’arguments graphiques.

par(mar = c(4,4,0,4))

# Line plot
tun <- world[world$Name == "Tunisia",]
alg <- world[world$Name == "Algeria",]
mar <- world[world$Name == "Morocco",]

plot(alg$Year, # Abscisses 
     alg$YOUNG_RT, # Ordonnées
     type = "l", # Type lignes
     ylim = c(20, 50), # Bornes min/max des ordonnées
     cex = .6, # Taille des points
     col = "blue", # Couleur de la ligne
     cex.lab = 0.7, # Réduit les labels d'un facteur de 0.7
     cex.axis = 0.6, # Réduit les labels des graduations d'un facteur de 0.6 
     xlab = "Années", # Label abscisses 
     ylab = "Part des jeunes (%) dans la population totale") # Label ordonnées

lines(mar$Year, 
      mar$YOUNG_RT, # Rajouter une ligne (Maroc)
      col = "darkgreen", 
      type = "l",
      cex = .6)

lines(tun$Year,  # Et Tunisie
      tun$YOUNG_RT, 
      type = "l", 
      cex = .6,
      col = "red")

# Organisation de la légende
legend(x = "topright",
       legend = c("Algérie", "Maroc", "Tunisie"),
       col = c("blue", "darkgreen", "red"), lty = 1,
       cex = .8)

Graphiques interactifs

Le package plotly (Sievert et al. 2024) permet d’intégrer une dimension interactive aux représentations graphiques. Pour utiliser les fonctions de ce package, il faut bien avoir en tête le format de données attendues (format long).

library(plotly)

# Importer jeu de données d'exemple de gapminder
df <- read.csv("data/world/gapminder.csv")

head(df)
      country continent year lifeExp      pop gdpPercap
1 Afghanistan      Asia 1952  28.801  8425333  779.4453
2 Afghanistan      Asia 1957  30.332  9240934  820.8530
3 Afghanistan      Asia 1962  31.997 10267083  853.1007
4 Afghanistan      Asia 1967  34.020 11537966  836.1971
5 Afghanistan      Asia 1972  36.088 13079460  739.9811
6 Afghanistan      Asia 1977  38.438 14880372  786.1134
# Créer un graphique de type "Multiple Trace Animations"
fig <- df %>%
  plot_ly(
    x = ~gdpPercap, 
    y = ~lifeExp, 
    size = ~pop, 
    color = ~continent, 
    frame = ~year, 
    text = ~country, 
    hoverinfo = "text",
    type = 'scatter',
    mode = 'markers'
  )

# Paramétrer l'échelle, les labels, etc.
fig <- fig %>% layout(
    xaxis = list(type = "log", title = "PIB par habitant"),
    yaxis = list(title = "Espérance de vie")
  )


fig
Pour aller plus loin

Consulter r-graph-gallery qui présente les possibilités graphiques les plus courantes, en langage de base R. L’exploration de la syntaxe ggplotet son package éponyme ggplot2 peut être utile pour la recherche de représentations graphiques avancées. Cette syntaxe est un peu différente du langage de base R et repose sur les principes de la “grammaire des graphiques”. Plusieurs manuels très bien construits permettent de rentrer dans l’univers de la visualisation de données avec ggplot2, comme Modern Data Visualization with R (Kabacoff 2023)

Les fournisseurs de données institutionnels distribuent parfois des tableaux de données peu adaptés à l’intégration dans un logiciel quel qu’il soit. La mise en forme des données dans un langage de programmation permet :

  • d’éviter d’avoir à manipuler ces fichiers manuellement, source d’erreurs.
  • d’avoir un façon générique et rapide d’importer un grand nombre de tableaux de données (et éviter de perdre de nombreuses heures de travail).

Voici un exemple adapté aux données de chômage par gouvernorats tunisiens produit par l’Institut National de la Statistique tunisien (INS). Le fichier brut ressemble à ceci :

Un court bloc de code permet de réorganiser le fichier comme désiré, et pourrait être utilement étendu aux autres fichiers distribués par l’INS.

# Import
df <- read_excel("data/tun/chomage.xls", # Chemin du fichier
                 sheet = "Sheet1", # Nom de la feuille Excel
                 skip = 3, # Retirer les trois premières lignes
                 col_types = c(rep("text", 3), rep("numeric", 6))) # Format des colonnes 

# Table de passage Recensement / codes internationaux
gouv <- read.csv("data/tun/table_GADM_TUN.csv")

# Reformatage
df <- data.frame(df)
df[,1:2] <- NULL # Retirer les deux premières colonnes

# Noms de colonnes
cols <- c("CHOM_INS_NON_", "CHOM_INS_PRI_", "CHOM_INS_SEC_", "CHOM_INS_SUP_",
          "CHOM_INS_UNK_", "CHOM_INS_TOT_") 
colnames(df)[1] <- "id_TUN"

# Chômeurs
# Hommes
df1 <- df[c(1:25),]
names(df1)[2:length(df1)] <- paste0(cols, "M")
gouv <- merge(gouv, df1, by = "id_TUN", all.x = TRUE)

# Femmes
df2 <- df[c(26:50),]
names(df2)[2:length(df2)] <- paste0(cols, "F")
gouv <- merge(gouv, df2, by = "id_TUN", all.x = TRUE)

# Total
df3 <- df[c(51:75),]
names(df3)[2:length(df3)] <- paste0(cols, "T")
gouv <- merge(gouv, df3, by = "id_TUN", all.x = TRUE)

head(gouv)
     id_TUN id_GADM CHOM_INS_NON_M CHOM_INS_PRI_M CHOM_INS_SEC_M CHOM_INS_SUP_M
1    ARIANA   TN.AN            484           2711           6876           2979
2      BEJA   TN.BJ           1211           2701           4577           1399
3 BEN AROUS   TN.BA            400           3496           8869           3512
4   BIZERTE   TN.BZ           1276           4592           8087           2568
5     GABES   TN.GB            305           2465           4890           2620
6     GAFSA   TN.GF            488           3031           7571           3576
  CHOM_INS_UNK_M CHOM_INS_TOT_M CHOM_INS_NON_F CHOM_INS_PRI_F CHOM_INS_SEC_F
1             69          13119            469           1459           4674
2             77           9965           1025           1319           2728
3             81          16358            696           2024           5960
4             74          16597            698           1775           4063
5            100          10380            340           1155           3007
6             67          14733            531           1711           4516
  CHOM_INS_SUP_F CHOM_INS_UNK_F CHOM_INS_TOT_F CHOM_INS_NON_T CHOM_INS_PRI_T
1           6610             44          13256            953           4170
2           2718             51           7841           2236           4020
3           6977             37          15694           1096           5520
4           5311             64          11911           1974           6367
5           7235            100          11837            645           3620
6           7639             80          14477           1019           4742
  CHOM_INS_SEC_T CHOM_INS_SUP_T CHOM_INS_UNK_T CHOM_INS_TOT_T
1          11550           9589            113          26375
2           7305           4117            128          17806
3          14829          10489            118          32052
4          12150           7879            138          28508
5           7897           9855            200          22217
6          12087          11215            147          29210

Une fois mis en forme, on peut réaliser un graphique sans avoir modifié le fichier initial provenant de l’INS.

# 2 graphiques par ligne
par(mfrow = c(1,2), mar = c(2,2,2,2))

# Boxplot hommes
boxplot(gouv$CHOM_INS_NON_M, # Chômeurs hommes par niveau d'instruction
        gouv$CHOM_INS_PRI_M, 
        gouv$CHOM_INS_SEC_M,
        gouv$CHOM_INS_SUP_M, 
        ylim = c(0, 12000), # Bornes min/max de l'axe des ordonnées
        main = "Hommes", # Titre plot
        names = c("Rien", "Primaire", "Secondaire", "Tertiaire"), # Labels (X)
        ylab = "Nombre de chômeurs, 2014", # Label (Y)
        col = "#adcaf7", # Couleur des box-plots
        cex.axis = .6, # Taille des labels des axes (réduit de 70 %)
        cex.title = .6) # Taille du label du titre (réduit de 70 %)

# Boxplot femmes
boxplot(gouv$CHOM_INS_NON_F, 
        gouv$CHOM_INS_PRI_F, 
        gouv$CHOM_INS_SEC_F,
        gouv$CHOM_INS_SUP_F,
        ylim = c(0, 12000),
        main = "Femmes",
        names = c("Rien", "Primaire", "Secondaire", "Tertiaire"), 
        col = "#ed9fb0", 
        cex.axis = .6, 
        cex.title = .6)

Nombre de chômeurs dans les gouvernorats tunisiens par niveau d’éducation

Export

write.csv() exporte un data.frame selon un chemin spécifié.

write.csv(x = gouv, # Objet à exporter
          dsn = "data/tun/gouv_chom.csv", # Chemin de fichier 
          row.names = FALSE) # Pour retirer les numéros de ligne


Données vectorielles

Le package sf permet d’importer et manipuler des couches géographiques vectorielles (points, lignes, polygones).

Import

Le geopackage est un format de données ouvert qui permet de stocker plusieurs couches géographiques dans un même fichier. La fonction st_layers() permet d’avoir un aperçu des couches présentes dans un fichier geopackage.

library(sf)

st_layers("data/tun/geom/tun_admin.gpkg")
Driver: GPKG 
Available layers:
   layer_name geometry_type features fields                      crs_name
1  delegation Multi Polygon      266     11 ETRS89-extended / LAEA Europe
2 gouvernorat Multi Polygon       24      4 ETRS89-extended / LAEA Europe
3      region Multi Polygon        6      2 ETRS89-extended / LAEA Europe

Importer des couches géographiques avec st_read(). Il s’agit ici des couches géographiques créées il y a quelques années par l’UAR RIATE.

del <- st_read("data/tun/geom/tun_admin.gpkg", layer = "delegation", quiet = TRUE)
gouv <- st_read("data/tun/geom/tun_admin.gpkg", layer = "gouvernorat", quiet = TRUE)
reg <- st_read("data/tun/geom/tun_admin.gpkg", layer = "region", quiet = TRUE)

Il s’agit bien d’objets sf.

class(del)
[1] "sf"         "data.frame"

La fonction st_read() peut aussi être employée pour des formats de fichiers .geojson, .shapefiles, etc.

reg <- st_read("data/tun/geom/map_reg.geojson", quiet = TRUE) 

Les objets sf sont des data.frame dont l’une des colonnes contient des géométries. Cette colonne est de la classe sfc (simple feature column) et chaque individu de la colonne est un sfg (simple feature geometry).

Ce format est très pratique dans la mesure où les données et les géométries sont intrinsèquement liées dans un même objet. Des informations additionnelles, propres aux couches d’informations géographiques sont aussi restituées (type de géométrie, emprise spatiale, système de coordonnées de référence).

head(del)
Simple feature collection with 6 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 4184398 ymin: 1282883 xmax: 4380989 ymax: 1539573
Projected CRS: ETRS89-extended / LAEA Europe
  del_code   del_nom_fr del_code_riate del_code_ins     del_nom_ar gou_code
1 TN.SF.AG       Agareb         TS2146        TN347          عقارب    TN.SF
2 TN.JE.AD   Ain Drahem         TS1224        TN225      عين دراهم    TN.JE
3 TN.SS.AK       Akouda         TS2115        TN316          اكودة    TN.SS
4 TN.KR.AL         Alaa         TS2216        TN417         العلاء    TN.KR
5 TN.BJ.AM       Amdoun         TS1212        TN213          عمدون    TN.BJ
6 TN.AN.AR Ariana Ville         TS1120        TN121 أريانة المدينة    TN.AN
   gou_nom gou_cap gou_cap_dist reg_code      reg_nom
1     Sfax       0         30.1       CE   Centre-est
2 Jendouba       0         38.5       NO   Nord-ouest
3   Sousse       0         14.3       CE   Centre-est
4 Kairouan       0         48.6       CO Centre-ouest
5     Beja       0         18.9       NO   Nord-ouest
6   Ariana       1          1.0       NE     Nord-est
                            geom
1 MULTIPOLYGON (((4375308 130...
2 MULTIPOLYGON (((4224769 152...
3 MULTIPOLYGON (((4373533 142...
4 MULTIPOLYGON (((4306097 141...
5 MULTIPOLYGON (((4243838 153...
6 MULTIPOLYGON (((4338370 153...

Aperçu des variables avec plot() :

plot(del)

Afficher juste les géométries:

# 3 cartes par ligne
par(mfrow = c(1,3), mar = c(2,2,2,2))

# Délégations
plot(del$geom, # Géométries uniquement
     col = "peachpuff", # Couleur du fond
     border = "white", # Couleur de bordure
     main = "Délégations") # Titre

# Gouvernorats
plot(gouv$geom, 
     col = "peachpuff", 
     border = "white",
     main = "Gouvernorats")

# "Régions"
plot(reg$geom, 
     col = "peachpuff",
     border = "white",
     main = "Régions")

Jointures attributaires

Nous pouvons joindre un data.frame à un objet sf en utilisant la fonction merge() et en s’appuyant sur des identifiants communs aux deux objets.

Attention à l’ordre des arguments, l’objet retourné sera du même type que x. Il n’est pas possible de faire une jointure attributaire en utilisant deux objets sf.

del <-  merge(
  x = del[,"del_code"],  # L'objet sf (seulement le champ del_code)
  y = del_df,          # le data.frame
  by.x = "del_code",  # identifiant dans x
  by.y = "del_code",  # identifiant dans y
  all.x = TRUE         # conserver toutes les lignes
)

head(del)
Simple feature collection with 6 features and 30 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 4321914 ymin: 1528813 xmax: 4346090 ymax: 1562156
Projected CRS: ETRS89-extended / LAEA Europe
  del_code         del_nom_fr     del_nom_ar gou_code gou_nom gou_cap
1 TN.AN.AR       Ariana Ville أريانة المدينة       AN  Ariana       1
2 TN.AN.ET    Cité Ettathamen     حي التضامن       AN  Ariana       0
3 TN.AN.KA kalaât El Andalous   قلعة الاندلس       AN  Ariana       0
4 TN.AN.LS             Soukra           سكرة       AN  Ariana       0
5 TN.AN.MN             Mnihla       المنيهلة       AN  Ariana       0
6 TN.AN.RA             Raoued           رواد       AN  Ariana       0
  gou_cap_dist reg_code  reg_nom popto_2004 popco_2004 immig_2004 emigr_2004
1          1.0       NE Nord-est      97687      97687      16961      15426
2          8.8       NE Nord-est      78311      78311       5651       5245
3         19.8       NE Nord-est      23045      15313        728        528
4          6.7       NE Nord-est      89151      89151      19129       2832
5          8.0       NE Nord-est      53752      40176       7258        985
6          8.8       NE Nord-est      60896      53911      14053       1443
  mobil_2004 menag_2004 ordin_2004 porta_2004 telef_2004 popto_2014 popco_2014
1      32387      27468       9751      22524      18596     114486     114486
2      10896      11950        430       4505       3824      84312      84312
3       1256       4709        188       1865        829      26796      18211
4      21961      21590       2979      13321       9262     129693     129693
5       8243      17119        685       7087       7019      89884      58641
6      15496      14276       2327       9022       6339     106414      94961
  immig_2014 emigr_2014 mobil2014 menag_2014 ordin_2014 porta_2014 telef_2014
1      15637      20448     36085      32498      25474      32308      19942
2       5028       6752     11780      22087       6836      21715       3496
3       1104        752      1856       6554       1701       6305        428
4      20786       5528     26314      33981      18191      33683       9831
5      14400       1828     16228      22781       8326      22531       4249
6      20128       2521     22649      27574      14284      27338       7505
  popto_2010 surfa_2010 idr_2011                       geometry
1     109500     18.612    0.638 MULTIPOLYGON (((4338370 153...
2      82922      3.376    0.386 MULTIPOLYGON (((4330561 153...
3      24367    188.206    0.383 MULTIPOLYGON (((4333853 156...
4     107829     27.895    0.557 MULTIPOLYGON (((4342368 153...
5      69247     22.907    0.466 MULTIPOLYGON (((4332840 153...
6      77419     57.080    0.531 MULTIPOLYGON (((4338691 154...
plot(del[,"idr_2011"])

Sélectionner des lignes, des colonnes

Les objets sf sont des data.frame, on peut donc sélectionner leur lignes et leur colonnes de la même manière.

# Sélection de lignes
sou <- del[del$gou_nom == "Sousse",]

# Sélection de colonnes
sou <- sou[,"idr_2011"]

# Ne conserver que les lignes avec une valeur
sou <- sou[!is.na(sou$idr_2011),]

plot(sou[,"idr_2011"])

Export d’une couche géographique

st_write("data/sousse_deleg.geojson")


Données raster

Le package terra (Hijmans 2024) permet d’importer et d’exporter des fichiers raster. Il repose sur la bibliothèque GDAL (GDAL/OGR contributors 2022) qui permet de lire et de traiter un très grand nombre de format d’images géographiques.

Import

La fonction rast() permet de créer et/ou d’importer des données raster. Nous importons ici un jeu de données au format .tif créé et distribué par WorldPop qui porte sur une estimation de la population tunisienne dans une résolution de 100 mètres.

library(terra)
pop <- rast("data/tun/raster/tun_ppp_2020_UNadj.tif")

Ce sont des objets de type SpatRaster.

pop
class       : SpatRaster 
dimensions  : 8788, 4881, 1  (nrow, ncol, nlyr)
resolution  : 0.0008333333, 0.0008333333  (x, y)
extent      : 7.530417, 11.59792, 30.23708, 37.56042  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : tun_ppp_2020_UNadj.tif 
name        : tun_ppp_2020_UNadj 

La fonction summary() est toujours utile pour un résumé statistique des cellules. Vu le nombre important de cellules, ce résumé est effectué sur un échantillon.

summary(pop)
 tun_ppp_2020_UNadj
 Min.   :  0.00    
 1st Qu.:  0.01    
 Median :  0.09    
 Mean   :  0.55    
 3rd Qu.:  0.25    
 Max.   :234.75    
 NA's   :49414     

Affichage

Par défaut et comme montré précédemment, la fonction plot() renvoie une légende continue sur les représentations cartographiques. On peut choisir de discrétiser cette information pour représenter des classes de valeur associée à une palette de couleurs personnalisée.

par(mfrow = c(1,2))

# Graphique en échelle continue
plot(pop)

# Avec discrétisation et paramétrage des couleurs
cuts <- c(0, 1, 2, 4, 8, 16, 32, 64, 235)
cols <- colorRampPalette(c("yellow", "darkgoldenrod1", "brown1"))
plot(pop, breaks = cuts, col = cols(8))

Modification de la zone d’étude

Le découpage d’un raster en fonction de l’étendue d’un autre objet, SpatVector ou SpatRaster, est réalisable avec la fonction crop(). Les deux couches de données doivent être dans la même projection.

# Transformer en WGS 84
sou <- st_transform(sou, 4326)

# Crop avec la délégation de Sousse
pop_sou <- crop(pop, sou)

# Plot
plot(pop_sou, breaks = cuts, col = cols(8))
plot(sou$geometry, col = NA, add = TRUE)

Export

La fonction writeRaster() permet d’enregistrer un objet SpatRaster.

writeRaster(x = pop_sou, filename = "data/tun/raster/pop_sou.tif")


Packages de données spatiales

De nombreux packages distribuent des données géographiques. Ils interfacent souvent des API qui permettent d’interroger des données mises à disposition sur le Web, directement avec R.

L’intérêt des API est de se connecter directement aux bases de données des fournisseurs de données, et de disposer des dernières mises à jour des données.

Au niveau international

Au niveau national

  • Brésil
    • geobr (Pereira et Goncalves 2024) : fournit un accès facile aux séries de données spatiales officielles du Brésil pour différentes années et découpage administratifs.
  • Chili -chilemapas (Vargas 2022) : donne accès aux divisions politiques et administratives du Chili.
  • Espagne
    • mapSpain (Hernangómez 2024b) : propose les limites administratives de l’Espagne à plusieurs niveaux (Communautés autonomes, Provinces, Municipalités), ainsi que des tuiles.
  • États-Unis
    • tidycensus (Walker et Herman 2024) : permet de charger des données et géométries du recensement américain en format sf et tidyverse
    • tigris (Walker 2024b) : donne accès aux éléments cartographiques fournis par le US Census Bureau TIGER, y compris les limites cartographiques, les routes et l’eau.
    • FedData (Bocinsky 2024) : automatise le téléchargement de données géospatiales disponibles à partir de plusieurs sources de données fédérées.
    • acs (Glenn 2019) : permet de télécharger et manipuler les données de l’American Community Survey et les données décennales du recensement des États-Unis.
    • censusapi (Recht 2022) : wrapper pour les API du Census Bureau des États-Unis.
    • idbr (Walker 2024a) : interface avec l’API de la base de données internationale du US Census Bureau.
    • ipumsr (Greg Freedman Ellis, Derek Burk, et Finn Roberts 2024) : Permet d’importer des données de recensement, d’enquête et géographiques fournies par l’IPUMS.
    • totalcensus (Li 2021) : permet d’extraire les données du recensement décennal et de l’American Community Survey au niveaux des block, block group et tract.
  • Finland
    • mapsFinland (Haukka 2023) : donne un accès à des cartes et données concernant la Finland.
  • France
  • Pologne
  • Uruguay
    • geouy (Detomasi 2023) : permet le chargement d’informations géographiques sur l’Uruguay.

geodata

Ce package facilite l’accès à des données géographique de référence sur le climat, la couverture du sol, les limites administratives et plusieurs autres jeux de données de référence au niveau mondial.

Function Description
bio_oracle Marine data from Bio-Oracle
cmip6_world, cmip6_tile Downscaled and calibrated CMIP6 projected future climate data
country_codes Country codes
crop_calendar_sacks Crop calendar data by Sacks et al
crop_monfreda Crop area and yield data for 175 crops by Monfreda et al.
crop_spam MapSPAM crop data (area, yield, value)
cropland Cropland density for the world derived from different sources (ESA, GLAD, QED)
elevation_3s, elevation_30s, elevation_global Elevation data
gadm, world Administrative boundaries for any country, or the entire world from GADM
landcover Landcover data derived from ESA WorldCover
footprint Human footprint data from the Last of the Wild project
osm OpenStreetMap data by country (places and roads)
population Population density Gridded Population of the World
soil_af Chemical and physical soil properties data for Africa for different soil depths
soil_af_elements Connect to or download chemical soil element concentration (for the 0-30 cm topsoil) data for Africa
soil_af_water Physical soil properties data for Africa for water balance computation
soil_af_isda Soil data for Africa derived from the iDSA data set
soil_world_vsi Virtually connect to the global SoilGrids data
soil_world Global soils data from SoilGrids
sp_occurrence Species occurrence records from the Global Biodiversity Information Facility
travel_time Travel time to and from cities and ports by Nelson et al.
worldclim_global, worldclim_country, worldclim_tile WorldClim glocal climate data

Extraction des couches géographiques d’altitude et de température en Tunisie.

library(geodata)
elev <- elevation_30s(country = "TUN", path = tempdir())
temp <- worldclim_country(country = "Tunisia", 
                          res = 10, 
                          var = "tavg",
                          path = tempdir())

Ce sont des objets de type SpatRaster.

class(temp)
[1] "SpatRaster"
attr(,"package")
[1] "terra"
class(elev)
[1] "SpatRaster"
attr(,"package")
[1] "terra"
# Les afficher
par(mfrow = c(1,2))

# Altitude
cols <- colorRampPalette(c("#31ad37", "#f5f752", "#fca330", "#9c5903"))
plot(elev, main = "Altitude", col = cols(50))
     
# Température
cols <- colorRampPalette(c("#4575b4", "#91bfdb", "#e0f3f8", 
                           "#fee090", "#fc8d59", "#d73027"))
plot(temp$TUN_wc2.1_30s_tavg_5, main = "Températures, Mai (1970-2000)",
     col = cols(50))

Natural Earth et Banque Mondiale.

Objectif

Nous cherchons ici à créer une carte des émissions de CO² en Afrique en utilisant les données de la banque mondiale et le fond de carte Natural Earth. Rien besoin de télécharger, mais il faut une connexion internet.

Plusieurs fonds de carte des pays du Monde à différents niveaux de généralisation cartographique sont disponibles avec la package rnaturalearth. Nous sélectionnons uniquement les pays africains.

library(rnaturalearth)

# Import des pays
country <- ne_countries(type = "countries", # pays
                        scale = "small",  # niveau de généralisation
                        returnclass = "sf") # objet retourné


# Conversion en projection Mercator
country <- st_transform(country, crs = "EPSG:3857")

# Si on s'intéresse à l'Afrique (modèle carto)
afr <- country[country$continent == "Africa",]

Le package wbstatspermet d’interroger l’API de la base de données de la Banque Mondiale. On peut faire une recherche pour trouver le nom des tables qui répondent à une requête par mots-clés avec la fonction wb_search().

library(wbstats)

wb_search("CO2 emissions")
# A tibble: 20 × 3
   indicator_id         indicator                                 indicator_desc
   <chr>                <chr>                                     <chr>         
 1 EN.ATM.CO2E.CP.KT    CO2 emissions from cement production (th… "Carbon dioxi…
 2 EN.ATM.CO2E.FF.KT    CO2 emissions from fossil-fuels, total (… "Fossil fuel …
 3 EN.ATM.CO2E.FF.ZS    CO2 emissions from fossil-fuels (% of to… "Fossil fuel …
 4 EN.ATM.CO2E.GF.KT    CO2 emissions from gaseous fuel consumpt… "Carbon dioxi…
 5 EN.ATM.CO2E.GF.ZS    CO2 emissions from gaseous fuel consumpt… "Carbon dioxi…
 6 EN.ATM.CO2E.GL.KT    CO2 emissions from gas flaring (thousand… "Carbon dioxi…
 7 EN.ATM.CO2E.KD.GD    CO2 emissions (kg per 2010 US$ of GDP)    "Carbon dioxi…
 8 EN.ATM.CO2E.KT       CO2 emissions (kt)                        "Carbon dioxi…
 9 EN.ATM.CO2E.LF.KT    CO2 emissions from liquid fuel consumpti… "Carbon dioxi…
10 EN.ATM.CO2E.LF.ZS    CO2 emissions from liquid fuel consumpti… "Carbon dioxi…
11 EN.ATM.CO2E.PC       CO2 emissions (metric tons per capita)    "Carbon dioxi…
12 EN.ATM.CO2E.PP.GD    CO2 emissions (kg per PPP $ of GDP)       "Carbon dioxi…
13 EN.ATM.CO2E.PP.GD.KD CO2 emissions (kg per 2017 PPP $ of GDP)  "Carbon dioxi…
14 EN.ATM.CO2E.SF.KT    CO2 emissions from solid fuel consumptio… "Carbon dioxi…
15 EN.ATM.CO2E.SF.ZS    CO2 emissions from solid fuel consumptio… "Carbon dioxi…
16 EN.CO2.BLDG.ZS       CO2 emissions from residential buildings… "CO2 emission…
17 EN.CO2.ETOT.ZS       CO2 emissions from electricity and heat … "CO2 emission…
18 EN.CO2.MANF.ZS       CO2 emissions from manufacturing industr… "CO2 emission…
19 EN.CO2.OTHX.ZS       CO2 emissions from other sectors, exclud… "CO2 emission…
20 EN.CO2.TRAN.ZS       CO2 emissions from transport (% of total… "CO2 emission…

La fonction wb_data() extrait les tables identifiées pour pour construire l’indicateur d’émissions de CO² par habitant. La série statistique est disponible de 1960 à 2023.

Avertissement

Chaque package reposant sur une API a son fonctionnement propre. Il est donc conseillé de consulter attentivement la documentation associée à ces packages (souvent précise) pour apprendre à les utiliser.

# Sélection des indicateurs
my_indicators = c("pop" = "SP.POP.TOTL",
                  "co2" = "EN.ATM.CO2E.KT")

# Interroger l'API
wb <- wb_data(my_indicators, return_wide = FALSE)

head(wb)
# A tibble: 6 × 11
  indicator_id indicator      iso2c iso3c country  date   value unit  obs_status
  <chr>        <chr>          <chr> <chr> <chr>   <dbl>   <dbl> <chr> <chr>     
1 SP.POP.TOTL  Population, t… AF    AFG   Afghan…  2023 NA      <NA>  <NA>      
2 SP.POP.TOTL  Population, t… AF    AFG   Afghan…  2022  4.11e7 <NA>  <NA>      
3 SP.POP.TOTL  Population, t… AF    AFG   Afghan…  2021  4.01e7 <NA>  <NA>      
4 SP.POP.TOTL  Population, t… AF    AFG   Afghan…  2020  3.90e7 <NA>  <NA>      
5 SP.POP.TOTL  Population, t… AF    AFG   Afghan…  2019  3.78e7 <NA>  <NA>      
6 SP.POP.TOTL  Population, t… AF    AFG   Afghan…  2018  3.67e7 <NA>  <NA>      
# ℹ 2 more variables: footnote <chr>, last_updated <date>
# Années disponibles
unique(wb$date)
 [1] 2023 2022 2021 2020 2019 2018 2017 2016 2015 2014 2013 2012 2011 2010 2009
[16] 2008 2007 2006 2005 2004 2003 2002 2001 2000 1999 1998 1997 1996 1995 1994
[31] 1993 1992 1991 1990 1989 1988 1987 1986 1985 1984 1983 1982 1981 1980 1979
[46] 1978 1977 1976 1975 1974 1973 1972 1971 1970 1969 1968 1967 1966 1965 1964
[61] 1963 1962 1961 1960

Nous gardons la dernière année pour laquelle les émissions de CO² sont disponibles 2020. La fonction reshape() est utilisée pour transformer le data.frame dans un format attendu (long vers large).

# Sélectionner la dernière année disponible et les colonnes utiles
wb <- wb[wb$date %in% 2020,]
wb <- wb[,c("iso3c", "indicator_id", "value")]
wb <- data.frame(wb)

# Formatage (format long + concaténation code indicateur + année ref
wb <- reshape(wb, idvar = "iso3c", timevar = "indicator_id", direction = "wide")

# Renommage des colonnes et création de l'indicateurs émissions par hab
names(wb)[2:3] <- c("co2_2020", "pop_2020")
wb$co2_hab_2020 <- wb$co2_2020 * 1000 / wb$pop_2020

La jointure attributaire avec le fond de carte précédemment importé avec rnaturalearth permet ainsi de constater que tout est en place pour procéder à une représentation cartographique du phénomène.

# Jointure attributaire
afr <- merge(afr[,c("adm0_a3", "name_fr")], wb, 
             by.x = "adm0_a3", by.y = "iso3c", all.x = TRUE) 

# Carte
plot(afr[,"co2_hab_2020"])

OpenStreetMap

OpenStreetMap (OSM) est un projet de cartographie participative qui a pour but de constituer une base de données géographiques libre à l’échelle mondiale. OSM vous permet de voir, modifier et utiliser des données géographiques dans le monde entier.

Plusieurs packages permettre d’extraire, interroger et visualiser des données issues d’OSM :

Objectif

Dans un premier temps nous cherchons à extraire d’OSM les géométries suivantes en Tunisie :

  • Découpages administratifs : Délégations et secteurs.
  • Équipements : objets décrits avec une clé OSM amenity (équipements utiles et importants).

A partir de ces données extraites, nous proposons ensuite de calculer les temps de parcours routiers des centres de secteur vers l’hôpital ou la clinique la plus accessible.

Extraction et préparation des données

Le package osmextract (Gilardi et Lovelace 2023) permet d’extraire des données depuis une base de données OSM directement et travailler en local sur des volumes de données très importants et ainsi d’éviter de surcharger un serveur Overpass turbo (utilisé par le package osmdata).

Découpages administratifs

La fonction oe_get() permet de télécharger un extrait de la base de données OSM pour une zone particulière et un type d’objet géographique. L’argument place correspond au nom du fichier *.pbf accessible sur le site Geofabrik. L’argument extra_tag permet de sélectionner les objets de la base OSM correspondant à une clé particulière (se référer à la documentation d’OSM pour choisir les clés).

On commence par extraire tous les polygones et les points OSM inclus en Tunisie. L’import d’extra_tags permet d’obtenir des informations supplémentaires qui serviront à filtrer la base de données.

library(osmextract)
osm_poly <- oe_get(place = "Tunisia",
                   layer = "multipolygons",
                   extra_tags = c("amenity", "ref:tn:hasc_2", "ref:tn:codegeo",
                                  "name:fr", "name:ar"))

osm_pt <- oe_get(place = "Tunisia",
                 layer = "points",
                 extra_tags = c("amenity", "name:fr", "name:ar"))

On peut ensuite filtrer l’extraction en fonction des objets / champs d’intérêt.

# Délégations et secteurs
admin <- osm_poly[!is.na(osm_poly$admin_level),] # Retirer pas d'attribut de niveau hiérarchique
deleg <- admin[admin$admin_level == 5,] # Délégations
sect <- admin[admin$admin_level == 6,] # Secteurs

# Ne garder que les champs utiles
deleg <- deleg[,c("osm_id", "ref_tn_codegeo", "ref_tn_hasc_2", "name_ar",  
                  "name_fr", "admin_level")]
sect <- sect[,c("osm_id", "ref_tn_codegeo", "name_ar",  
                  "name_fr", "admin_level")]

Ce bloc de code, dont on ne détaillera pas le contenu ici, permet d’harmoniser le nom des champs entre les couches géographiques et de disposer des délégations et gouvernorats d’appartenance des couches géographiques des secteurs et délégations.

# Délégation d'appartenance du secteur
sect_pt <- st_centroid(sect) # Centroid du secteur
sect_pt <- st_intersection(sect_pt, deleg) # Intersection avec couche délégation
sect_pt <- st_set_geometry(sect_pt, NULL) # Retirer géométries
sect <- merge(sect, # Enrichir les secteurs du code d'appartenance de la délégation
              sect_pt[,c("osm_id", "ref_tn_hasc_2")], 
              by = "osm_id", 
              all.x = TRUE)
sect <- sect[!duplicated(sect$osm_id),] # Retirer valeurs dupliquées

# Renommer colonnes
names(sect)[2] <- "id_tn"
names(sect)[6] <- "id_hasc_deleg"
sect$id_hasc_gouv <- substr(sect$id_hasc_deleg, 1, 5) # Gouvernorat d'appartenance
names(sect)[2,6] <- (c("tn_codegeo", "int_codegeo"))
names(deleg)[2:3] <- c("id_tn", "id_hasc_deleg")
deleg$id_hasc_gouv <- substr(deleg$id_hasc_deleg, 1, 5)

Équipements

Les objets décrits par la clé amenity sont de nature hétérogènes (points, polygones). Ce bloc de code extrait les centroides des polygones et les associent à la couche de points initiale.

Ces couches géographiques sont finalement exportées dans un geopackage tun_osm.gpkg.

# Consolidation des géométries
sel_poly <- osm_poly[!is.na(osm_poly$amenity),] # Retirer les objets qui n'ont pas de clé amenity
sel_pt <- osm_pt[!is.na(osm_pt$amenity),]
sel_poly <- st_make_valid(sel_poly) # Consolider les géométries des polygones
sel_pt2 <- st_centroid(sel_poly) # extraire le centroide
cols <- intersect(names(sel_pt), names(sel_pt2)) # Garder les colonnes identiques
sel_pt <- rbind(sel_pt[,cols, sel_pt2[,cols]]) # Combiner points et polygones
sel_pt <- sel_pt[,c("osm_id", "amenity", "name_ar", "name_fr")]

# Exporter les couche ainsi créées
st_write(sel_pt, "data/tun/geom/tun_osm.gpkg", layer = "poi")
st_write(deleg, "data/tun/geom/tun_osm.gpkg", layer = "deleg")
st_write(sect, "data/tun/geom/tun_osm.gpkg", layer = "sect")
Pour aller plus loin

Ce document (Giraud 2017) montre comment créer un fond de carte avec des données OSM, extraire des objets d’intérêt (bars et restaurants) avec le package osmdata. Il propose également plusieurs pistes cartographiques pour manipuler et visualiser ces données. Notons simplement que la cartographie utilise les fonctions du package cartography, plus maintenu, et qu’il est conseillé d’utiliser dorénavant le package mapsf.

Résultats de l’extraction

Dénombrement des 15 objets OSM les plus fréquents décrits par la clé amenity en Tunisie.

# Import des aménités préparées en amont
osm_pt <- st_read("data/tun/geom/tun_osm.gpkg", layer = "poi", quiet = TRUE)

# Nombre de points avec le tag "amenity". 
tbl <- table(osm_pt$amenity)
tbl <- tbl[order(tbl, decreasing = TRUE)]
tbl[1:15]

            cafe       restaurant place_of_worship             bank 
            2680             1272             1094              996 
          school         pharmacy        fast_food             fuel 
             916              700              597              525 
     post_office           police              bar          parking 
             327              153              135              135 
         doctors              atm          dentist 
             113              109              109 
# On ne prend ici que les cliniques et hôpitaux
osm_pt <- osm_pt[osm_pt$amenity %in% c("clinic", "hospital"), ]

Nombre de délégations, secteurs et cliniques/hôpitaux inclus dans OSM.

# Import des unités territoriales préparées en amont
deleg <- st_read("data/tun/geom/tun_osm.gpkg", layer = "deleg", quiet = TRUE)
sect <- st_read("data/tun/geom/tun_osm.gpkg", layer = "sect", quiet = TRUE)

# Nombre d'objets (lignes)
nrow(deleg)
[1] 266
nrow(sect)
[1] 2104
nrow(osm_pt)
[1] 119

Le résultat sous forme d’une carte interactive avec le package leaflet, qui sera présenté ultérieurement dans la formation.

Code
# Choisir la palette
library(leaflet)
pal <- colorFactor(palette = c("red", "gold"),
                   domain = c("clinic", "hospital"))

# Cartographie interactive
leaflet(osm_pt) |> # Emprise = hôpitaux
  addProviderTiles("OpenStreetMap.HOT") |> # Type de tuiles chargées
  addPolygons(data = sect, # Secteurs
              col = "white",
              fillColor = "lightgrey",
              fillOpacity = 0.7,
              weight = 1,
              popup = paste0("<b>", sect$id_tn, "<br></b>",
                             sect$name_ar, "<br>", sect$name_fr),
              group = "Secteurs")|>
  addPolygons(data = deleg, # Délégations
              col = "darkgrey",
              fill = "lightgrey",
              fillOpacity = 0,
              weight = 1.2,
              popup = paste0("<b>", deleg$id_tn, " / ",
                             deleg$id_hasc_deleg,"<br></b>",
                             deleg$name_ar, "<br>", deleg$name_fr),
              group = "Délégations") |>
  addCircleMarkers(radius = 4, # Hôpitaux
                   stroke = FALSE,
                   color = ~ pal(amenity),
                   fillOpacity = 1,
                   popup = paste0(osm_pt$name_ar, "<br>", osm_pt$name_fr),
                   group = "Cliniques et hôpitaux") |>
  addLegend(pal = pal, # Légende pour différencier cliniques et hôpitaux
            values = c("clinic", "hospital"),
            opacity = 0.7,
            title = "OSM amenity",
            position = "bottomright") |>
  addLayersControl(overlayGroups = c("Secteurs", "Délégations", "Cliniques et hôpitaux"),
                   options = layersControlOptions(collapsed = FALSE)) 

Matrices de temps et itinéraires

Le package osrm (Giraud (2022)) sert d’interface entre R et le service de calcul d’itinéraire OSRM (Luxen et Vetter, 2011). Ce package permet de calculer des matrices de temps et de distances, des itinéraires routiers, des isochrones. Le package utilise par défaut le serveur de démo d’OSRM. En cas d’utilisation intensive il est fortement recommandé d’utiliser sa propre instance d’OSRM avec Docker.

La fonction osrmTable() permet de calculer des matrices de distances ou de temps par la route. Nous effectuons cette opération entre le centroide de chaque secteur et les cliniques / hôpitaux extraits plus haut avec osmextract.

# Choisir un gouvernorat
sel <- sect[sect$id_hasc_gouv == "TN.SS",]
sel <- sel[!is.na(sel$osm_id),]
sel <- sel[!duplicated(sel$osm_id),]

# Origines (centroides des secteurs)
ori <- st_centroid(sel)

# Considérer les cliniques et hôpitaux dans un voisinage de 20 km autour de la délégation
osm_pt <- st_transform(osm_pt, 2088) # Transformer en coordonnées planaires
sel <- st_transform(sel, 2088)
dest <- st_filter(osm_pt, st_buffer(sel, 20000), .predicate = st_intersects)
dest <- st_transform(dest, 4326)

# Calcul de temps de trajets avec OSRM (pas grosse requête)
library(osrm)
df <- osrmTable(src = ori, dst = dest, measure = "duration")

# Extraire les temps de trajet
df <- data.frame(df$durations)

# Formater la table d'une manière arrangeante
colnames(df) <- as.character(dest$osm_id)
row.names(df) <- as.character(ori$osm_id)
head(df)
        2109034339 2981447790 3756156059 4543660394 4967820121 4971604621
7105249       55.4       68.3       58.6       58.7       80.3       78.4
7105250       58.2       71.1       61.4       61.5       83.1       81.2
7105251       56.9       69.8       60.1       60.2       81.8       79.9
7105252       54.0       66.9       57.2       57.3       78.9       77.0
7105253       45.1       58.0       48.3       48.4       70.0       68.1
7105254       58.7       71.6       62.0       62.1       83.6       81.7
        5516163723 5673214937 5778893454 7919870954 10759343505 11756738948
7105249       59.5       89.7       24.6       79.2        60.6        60.1
7105250       62.3       92.5       25.8       82.5        63.4        62.9
7105251       61.0       91.2       27.4       83.0        62.2        61.6
7105252       58.1       88.3       33.2       77.0        59.3        58.7
7105253       49.2       79.4       29.9       68.1        50.3        49.8
7105254       62.9       93.1       48.4       77.6        64.0        63.4

On peut ensuite assez aisément avec la fonction apply() extraire le temps de trajet minimal par secteur.

time <- apply(df, 1, min) 
time <- data.frame(time)
time$osm_id <- row.names(time)
head(time)
        time  osm_id
7105249 24.6 7105249
7105250 25.8 7105250
7105251 27.4 7105251
7105252 33.2 7105252
7105253 29.9 7105253
7105254 48.4 7105254

Nous finissons par une jointure attributaire avec la couche géographique des secteurs, ce qui permet de réaliser une représentation cartographique des résultats avec le package mapsf qui sera présenté ultérieurement dans la formation.

Code
sel <- merge(sel, time, by = "osm_id", all.x = TRUE)
deleg <- st_transform(deleg, 2088)
dest <- st_transform(dest, 2088)

library(mapsf)
par(mfrow = c(1,1))
mf_init(sel)
mf_map(deleg, col = "lightgrey", 
       border = NA,
       add = TRUE)
mf_map(sel, 
       type = "choro",
       var = "time",
       nbreaks = 4,
       border = "white", 
       leg_title = "Minutes en voiture", 
       add = TRUE)
mf_map(deleg, col = NA, 
       border = "black",
       add = TRUE)
mf_map(dest, pch = 21, col = NA, bg = "red", add = TRUE)
mf_scale(size = 10)
mf_title("Délégation de Sousse : Temps de trajet vers l'hôpital ou la clinique la plus proche")
mf_credits(paste0("Source: © OpenStreetMap et Contributeurs, 2024\n",
                  "NB/ Toute chose égale par rapport à la complétude d'OSM. Contribuez pour compléter la carte le cas échéant !"))

A vous de jouer !

L’ensemble des données présentées dans ce tutoriel ont été rassemblées dans le dossier data. Vous êtes invités sur un espace de votre choix (Monde / Tunisie) à :

  • Importer un jeu de données tabulaire.
  • Importer une couche géographique vectorielle.
  • Faire une jointure attributaire.
  • Choisir une variable et réaliser un résumé statistique (summary()) et un graphique présentant la distribution des valeurs de la variable.
  • Faire une carte avec la fonction plot()

Session Info

sessionInfo()
R version 4.4.0 (2024-04-24 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)

Matrix products: default


locale:
[1] LC_COLLATE=French_France.utf8  LC_CTYPE=French_France.utf8   
[3] LC_MONETARY=French_France.utf8 LC_NUMERIC=C                  
[5] LC_TIME=French_France.utf8    

time zone: Europe/Paris
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] mapsf_0.9.0         osrm_4.1.1          wbstats_1.0.4      
 [4] rnaturalearth_1.0.1 terra_1.7-71        plotly_4.10.4      
 [7] ggplot2_3.5.1       readxl_1.4.3        leaflet_2.2.2      
[10] sf_1.0-16          

loaded via a namespace (and not attached):
 [1] gtable_0.3.5            xfun_0.43               htmlwidgets_1.6.4      
 [4] tzdb_0.4.0              leaflet.providers_2.0.0 vctrs_0.6.5            
 [7] tools_4.4.0             crosstalk_1.2.1         generics_0.1.3         
[10] curl_5.2.1              tibble_3.2.1            proxy_0.4-27           
[13] fansi_1.0.6             pkgconfig_2.0.3         KernSmooth_2.23-22     
[16] data.table_1.15.4       RColorBrewer_1.1-3      lifecycle_1.0.4        
[19] compiler_4.4.0          farver_2.1.1            stringr_1.5.1          
[22] munsell_0.5.1           codetools_0.2-20        maplegend_0.1.0        
[25] htmltools_0.5.8.1       class_7.3-22            yaml_2.3.8             
[28] lazyeval_0.2.2          pillar_1.9.0            jquerylib_0.1.4        
[31] tidyr_1.3.1             classInt_0.4-10         wk_0.9.1               
[34] tidyselect_1.2.1        digest_0.6.35           stringi_1.8.3          
[37] dplyr_1.1.4             purrr_1.0.2             fastmap_1.1.1          
[40] grid_4.4.0              colorspace_2.1-0        cli_3.6.2              
[43] magrittr_2.0.3          utf8_1.2.4              e1071_1.7-14           
[46] RcppSimdJson_0.1.11     readr_2.1.5             withr_3.0.0            
[49] scales_1.3.0            rmarkdown_2.26          httr_1.4.7             
[52] cellranger_1.1.0        hms_1.1.3               evaluate_0.23          
[55] knitr_1.46              viridisLite_0.4.2       s2_1.1.6               
[58] rlang_1.1.3             isoband_0.2.7           Rcpp_1.0.12            
[61] glue_1.7.0              DBI_1.2.2               mapiso_0.3.0           
[64] rstudioapi_0.16.0       jsonlite_1.8.8          R6_2.5.1               
[67] units_0.8-5            

Les références

Arel-Bundock, Vincent. 2022. WDI: World Development Indicators and Other World Bank Data. https://CRAN.R-project.org/package=WDI.
Aybar, Cesar. 2023. rgee: R Bindings for Calling the ’Earth Engine’ API. https://CRAN.R-project.org/package=rgee.
Bocinsky, R. Kyle. 2024. FedData: Functions to Automate Downloading Geospatial Data Available from Several Federated Data Sources. https://CRAN.R-project.org/package=FedData.
Busetto, Lorenzo, et Luigi Ranghetti. 2016. « MODIStsp: an R package for preprocessing of MODIS Land Products time series ». Computers & Geosciences 97: 40‑48. https://doi.org/10.1016/j.cageo.2016.08.020.
Cambon, Jesse, Diego Hernangómez, Christopher Belanger, et Daniel Possenriede. 2021. tidycensus: Load US Census Boundary and Attribute Data as ’tidyverse’ and ’sf’-Ready Data Frames. https://CRAN.R-project.org/package=tidycensus.
Carteron, Paul. 2023. happign: R Interface to ’IGN’ Web Services. https://CRAN.R-project.org/package=happign.
Detomasi, Richard. 2023. « geouy: Geographic Information of Uruguay ». https://github.com/RichDeto/geouy.
Dyba, Krzysztof, et Jakub Nowosad. 2021. « rgugik: Search and Retrieve Spatial Data from the Polish Head Office of Geodesy and Cartography in R ». Journal of Open Source Software 6 (59): 2948. https://doi.org/10.21105/joss.02948.
GDAL/OGR contributors. 2022. GDAL/OGR Geospatial Data Abstraction software Library. Open Source Geospatial Foundation. https://doi.org/10.5281/zenodo.5884351.
Gilardi, Andrea, et Robin Lovelace. 2023. osmextract: Download and Import Open Street Map Data Extracts. https://CRAN.R-project.org/package=osmextract.
Giraud, Timothée. 2017. Cartographic Explorations of the OpenStreetMap Database with R. https://rcarto.github.io/caRtosm/index.html.
———. 2022. « osrm: Interface Between R and the OpenStreetMap-Based Routing Service OSRM ». Journal of Open Source Software 7 (78): 4574. https://doi.org/10.21105/joss.04574.
———. 2024. maptiles: Download and Display Map Tiles. https://CRAN.R-project.org/package=maptiles.
Giraud, Timothée, et Hugues Pecout. 2024. Géomatique avec R. https://doi.org/https://doi.org/10.5281/zenodo.5906212.
Glenn, Ezra Haber. 2019. acs: Download, Manipulate, and Present American Community Survey and Decennial Data from the US Census. https://CRAN.R-project.org/package=acs.
Greg Freedman Ellis, Derek Burk, et Finn Roberts. 2024. ipumsr: An R Interface for Downloading, Reading, and Handling IPUMS Data. https://CRAN.R-project.org/package=ipumsr.
Haukka, Jari. 2023. mapsFinland: Maps of Finland. https://CRAN.R-project.org/package=mapsFinland.
Hernangómez, Diego. 2024a. giscoR: Download Map Data from GISCO API - Eurostat (version 0.4.2). https://doi.org/10.5281/zenodo.4317946.
———. 2024b. mapSpain: Administrative Boundaries of Spain (version 0.9.0). https://doi.org/10.5281/zenodo.5366622.
Hernangómez, Diego, et Jindra Lacko. 2024. nominatimlite: Interface with ’Nominatim’ API Service. https://cran.r-project.org/web/packages/nominatimlite/index.html.
Hijmans, Robert J. 2024. terra: Spatial Data Analysis. https://CRAN.R-project.org/package=terra.
Hijmans, Robert J., Márcia Barbosa, Aniruddha Ghosh, et Alex Mandel. 2023. geodata: Download Geographic Data. https://CRAN.R-project.org/package=geodata.
Hollister, Jeffrey, Tarak Shah, Jakub Nowosad, Alec L. Robitaille, Marcus W. Beck, et Mike Johnson. 2023. elevatr: Access Elevation Data from Various APIs. https://doi.org/10.5281/zenodo.8335450.
Kabacoff, Robert. 2023. Modern Data Visualization with R. https://rkabacoff.github.io/datavis/.
Lahti, Leo, Janne Huovari, Markus Kainu, et Przemyslaw Biecek. 2017. « Retrieval and Analysis of Eurostat Open Data with the eurostat Package ». The R Journal 9 (1): 385‑92. https://doi.org/10.32614/RJ-2017-019.
Leclerc, Hadrien. 2022. insee: Tools to Easily Download Data from INSEE BDM Database. https://CRAN.R-project.org/package=insee.
Li, Guanglai. 2021. totalcensus: Extract Decennial Census and American Community Survey Data. https://CRAN.R-project.org/package=totalcensus.
Mark Padgham, Bob Rudis, Robin Lovelace, et Maëlle Salmon. 2017. « osmdata ». Journal of Open Source Software 2 (14): 305. https://doi.org/10.21105/joss.00305.
Massicotte, Philippe, et Andy South. 2023. rnaturalearth: World Map Data from Natural Earth. https://CRAN.R-project.org/package=rnaturalearth.
Pereira, Rafael H. M., et Caio Nogueira Goncalves. 2024. geobr: Download Official Spatial Data Sets of Brazil. https://CRAN.R-project.org/package=geobr.
Piburn, Jesse. 2020. wbstats: Programmatic Access to the World Bank API. Oak Ridge, Tennessee: Oak Ridge National Laboratory. https://doi.org/10.11578/dc.20171025.1827.
Ranghetti, Luigi, Mirco Boschetti, Francesco Nutini, et Lorenzo Busetto. 2020. « sen2r: An R toolbox for automatically downloading and preprocessing Sentinel-2 satellite data ». Computers & Geosciences 139: 104473. https://doi.org/10.1016/j.cageo.2020.104473.
Read, Jordan S., Jordan I. Walker, Alison Appling, David L. Blodgett, Emily K. Read, et Luke A. Winslow. 2015. « geoknife: Reproducible web-processing of large gridded datasets ». Ecography. https://doi.org/10.1111/ecog.01880.
Recht, Hannah. 2022. censusapi: Retrieve Data from the Census APIs. https://CRAN.R-project.org/package=censusapi.
Rowlingson, Barry. 2019. geonames: Interface to the "Geonames" Spatial Query Web Service. https://CRAN.R-project.org/package=geonames.
Sievert, Carson, Chris Parmer, Toby Hocking, Scott Chamberlain, et Karthik Ram. 2024. plotly: Create Interactive Web Graphics via ’plotly.js’. https://cran.r-project.org/web/packages/plotly/index.html.
Sparks, Adam H. 2018. « nasapower: A NASA POWER Global Meteorology, Surface Solar Energy and Climatology Data Client for R ». The Journal of Open Source Software 3 (30): 1035. https://doi.org/10.21105/joss.01035.
Vargas, Mauricio. 2022. chilemapas: Mapas de las Divisiones Politicas y Administrativas de Chile (Maps of the Political and Administrative Divisions of Chile). https://CRAN.R-project.org/package=chilemapas.
Walker, Kyle. 2024a. idbr: R Interface to the US Census Bureau International Data Base API. https://CRAN.R-project.org/package=idbr.
———. 2024b. tigris: Load Census TIGER/Line Shapefiles. https://CRAN.R-project.org/package=tigris.
Walker, Kyle, et Matt Herman. 2024. tidycensus: Load US Census Boundary and Attribute Data as ’tidyverse’ and ’sf’-Ready Data Frames. https://CRAN.R-project.org/package=tidycensus.
Watson, Oliver J, Rich FitzJohn, et Jeffrey W Eaton. 2019. « rdhs: an R package to interact with The Demographic and Health Surveys (DHS) Program datasets ». Wellcome Open Research 4: 103. https://doi.org/10.12688/wellcomeopenres.15311.1.
Weidmann, Nils B., Guy Schvitz, et Luc Girardin. 2021. cshapes: The CShapes 2.0 Dataset and Utilities. https://CRAN.R-project.org/package=cshapes.